//      cds_main.cc
//      
//      Copyright 2012 tqlong <tqlong@espada>
//      
//      This program is free software; you can redistribute it and/or modify
//      it under the terms of the GNU General Public License as published by
//      the Free Software Foundation; either version 2 of the License, or
//      (at your option) any later version.
//      
//      This program is distributed in the hope that it will be useful,
//      but WITHOUT ANY WARRANTY; without even the implied warranty of
//      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//      GNU General Public License for more details.
//      
//      You should have received a copy of the GNU General Public License
//      along with this program; if not, write to the Free Software
//      Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
//      MA 02110-1301, USA.

/* example run

../cds --fom --input=../cds.fasta.train --ll=../cds.fasta.test.fix.60 --ll=../intron.fasta.test.fix.60 --llout=pospos --llout=negpos --fom_stats=fom_stats --period=3 --order=1 --order=2 --order=3 --order=4 --order=5 --order=6 --order=7 --order=8 --imm_xi_ll --imm_T=100 --fom_ll
 */

//#include <iostream>
#include <boost/program_options.hpp>
#include <boost/date_time/posix_time/posix_time.hpp> 
#include <boost/foreach.hpp>

#include "helper/stl_helper.h"
#include "helper/ublas_helper.h"
#include "helper/ml_data.h"
#include "helper/rand_helper.h"
#include "helper/dna_processing.h"
#include "helper/suffix_array.h"

#include "cds_vom.h"

namespace po = boost::program_options;
po::options_description desc("Allowed options");
po::variables_map vm;

void process_options(int argc, char** argv);
void usage(int exitcode = 0);

void read_data();
void build_fom();
void build_new_fom();
void build_vom();
void fix_length();
void split_train_test();
void accuracy();

boost::posix_time::ptime time_start_all(boost::posix_time::microsec_clock::local_time());
boost::posix_time::ptime time_start(boost::posix_time::microsec_clock::local_time());

boost::posix_time::ptime now() 
{
  return boost::posix_time::microsec_clock::local_time();
}

boost::posix_time::time_duration duration() 
{
  boost::posix_time::time_duration d = now() - time_start;
  time_start = now();
  return d;
}

boost::posix_time::time_duration duration_all() 
{
  return now() - time_start_all;
}

int main(int argc, char** argv)
{
  //MSG("BOOST SUB-VERSION " << BOOST_VERSION / 100 % 100);

  process_options(argc, argv);

  if (vm.count("help"))
    usage();
  else if (vm.count("read"))
    read_data();
  else if (vm.count("fom"))
    build_fom();
  else if (vm.count("new_fom"))
    build_new_fom();
  else if (vm.count("vom"))
    build_vom();
  else if (vm.count("fix_length"))
    fix_length();
  else if (vm.count("train_test"))
    split_train_test();
  else if (vm.count("accuracy"))
    accuracy();
  else
    usage();

  MSG("Duration = " << duration_all());
}

void process_options(int argc, char** argv)
{
  desc.add_options()
    ("help", "produce this usage")
    ("input", po::value<std::string>(), "input file")
    ("validation", po::value<std::string>(), "validation file (DI)")
    ("output", po::value<std::string>(), "output file")
    ("verbose", "print debug information")

    // likelihood computation
    ("ll", po::value<std::vector<std::string> >(), "sequence files to compute log-likelihood")
    ("llout", po::value<std::vector<std::string> >(), "likelihood output files")
    ("fom_ll",  "Compute likelihood of FOM model")
    ("imm_xi_ll", "Compute likelihood of IMM model (chi-square)")
    ("imm_chi_ll", "Compute likelihood of IMM model (chi-square)")
    ("imm_prune_chi_ll", "Compute likelihood of IMM model (chi-squared then prune)")
    ("imm_prune_di_ll", "Compute likelihood of IMM model (deleted interpolation then prune)")
    ("prune_ll", "Compute likelihood of pruned model")
    ("prune_xi_ll", "Compute likelihood of pruned IMM model")
    ("imm_di_ll", "Compute likelihood of IMM model (deleted interpolation)")
    
    // Accuracy calculation
    ("negneg", po::value<std::string>(), "negative set, negative model output")
    ("negpos", po::value<std::string>(), "negative set, positive model output")
    ("posneg", po::value<std::string>(), "positive set, negative model output")
    ("pospos", po::value<std::string>(), "positive set, positive model output")

    // sequence set manipulation
    ("length", po::value<int>(), "fixed length to cut")
    ("ratio", po::value<double>()->default_value(0.5), "ratio of train set / all")

    // VOM param
    ("order", po::value<std::vector<int> >(), "FOM orders")
    ("Tmin", po::value<std::vector<int> >(), "VOM min threshold")
    ("period", po::value<int>(), "periodic (1, 3)")
    ("eta", po::value<std::vector<double> >(), "VOM pruned threshold")

    // IMM param
    ("imm_T", po::value<int>(), "Threshold for number of occurence (IMM chi-square)")
    ("imm_Ts", po::value<vi_t>(), "Threshold for number of occurence (IMM chi-squared)")
    ("imm_delta", po::value<std::vector<double> >(), "Threshold for number of occurence (IMM chi-squared)")

    // output options
    ("fom_stats", po::value<std::string>(), "fom statistics file")

    //JOB TYPE
    ("read", "read data file\n--input [--output]")
    ("fom", "Build FOM\n--input [--output] [--order]+ --periodic")
    ("vom", "Build VOM\n--input [--output] [--order]+ --periodic")
    ("fix_length", "split long sequence to fixed length\n--input [--output] --length")
    ("train_test", "split set of sequences to train set and test seth\n --input [--train] [--test] --ratio")
    ("accuracy", "--negneg --negpos --pospos --posneg --output")

    ("new_fom", "Build FOM\n--input [--output] [--order]+ --periodic")
    ;
  po::store(po::parse_command_line(argc, argv, desc), vm);
  po::notify(vm);
}

void usage(int exitcode)
{
  std::cout << desc << std::endl;
  exit(exitcode);
}

void read_data()
{
  std::string inFile = VM(vm, "input", std::string);

  Sequences S;
  S.fasta_read(inFile);
  MSG("read " << sz(S) << " sequences.");

  sastring_vec_t SA(sz(S));
  std::copy(all_element(S), SA.begin());
  MSG("suffix array of " << sz(SA) << " sequences built.");
}

void fix_length()
{
  std::string input = VM(vm, "input", std::string);
  int length = VM(vm, "length", int);
  std::string output = input + ".fix." + tos(length);
  if (vm.count("output")) output = VM(vm, "output", std::string);
  MSG("++ Split to short length input=" << input 
      << " output=" << output << " length=" << length);

  Sequences seqs;
  seqs.fasta_read(input);

  FileAccessWrite(out, output);
  seqs.split(out, length);
}

void split_train_test() 
{
  std::string input = VM(vm, "input", std::string);
  std::string train = input + ".train";
  std::string test = input + ".test";
  if (vm.count("train")) train = VM(vm, "train", std::string);
  if (vm.count("test")) test = VM(vm, "test", std::string);
  double ratio = VM(vm, "ratio", double);

  MSG("++ Split train test input=" << input << " ratio=" << ratio
      << " train=" << train << " test=" << test);

  Sequences seqs, trainSeqs, testSeqs;
  seqs.fasta_read(input);
  MSG("n_total=" << sz(seqs));

  RAND::randu_t randu = RAND::randu();
  for (int i = 0; i < sz(seqs); i++) {
    double r = randu();
    if (r < ratio) { // to train set
      trainSeqs.pb(seqs[i]);
      trainSeqs.annotation_.pb(seqs.annotation_[i]);
    }
    else { // to test set
      testSeqs.pb(seqs[i]);
      testSeqs.annotation_.pb(seqs.annotation_[i]);
    }
  }
  MSG("n_train=" << sz(trainSeqs));
  MSG("n_test=" << sz(testSeqs));

  FileAccessWrite(trainout, train);
  FileAccessWrite(testout, test);
  trainSeqs.print(trainout, true);
  testSeqs.print(testout, true);
}

void accuracy()
{
  std::string posposFile = VM(vm, "pospos", std::string);
  std::string negposFile = VM(vm, "negpos", std::string);
  std::string posnegFile = VM(vm, "posneg", std::string);
  std::string negnegFile = VM(vm, "negneg", std::string);
  std::string output = VM(vm, "output", std::string);
  MSG("++ Accuracy pospos=" << posposFile << " posneg=" << posnegFile
      << " negneg=" << negnegFile << " negpos=" << negposFile);

  FileAccessRead(pospos, posposFile);
  FileAccessRead(negpos, negposFile);
  FileAccessRead(posneg, posnegFile);
  FileAccessRead(negneg, negnegFile);

  // positive set
  int line = 1;
  double n = 0, n_pos = 0, n_correct = 0, n_pos_correct = 0;

  while (!pospos.eof()) {
    FileAccessLine(pospos, posposAnno);
    FileAccessLine(posneg, posnegAnno);
    //MSG("pp=" << posposAnno << " pn=" << posnegAnno);
    if (posposAnno.empty() || posnegAnno.empty()) {
      MSG("positive set detect end of file line=" << line);
      break;
    }
    if (posposAnno[0] != '>' || posnegAnno[0] != '>') 
      FATAL_ERROR("annotation not found line=" << line);
    line++;

    FileAccessLine(pospos, posposLLStr);
    FileAccessLine(posneg, posnegLLStr);
    
    double posposLL = s2<double>(posposLLStr);
    double posnegLL = s2<double>(posnegLLStr);

    n++;
    n_pos++;
    if (posposLL > posnegLL) {
      n_correct++;
      n_pos_correct++;
    }
    line++;
  }

  line = 1;
  while (!negneg.eof()) {
    FileAccessLine(negpos, negposAnno);
    FileAccessLine(negneg, negnegAnno);
    if (negposAnno.empty() || negnegAnno.empty()) {
      MSG("negative set detect end of file line=" << line);
      break;
    }
    if (negposAnno[0] != '>' || negnegAnno[0] != '>') 
      FATAL_ERROR("annotation not found line=" << line);
    line++;

    FileAccessLine(negpos, negposLLStr);
    FileAccessLine(negneg, negnegLLStr);
    
    double negposLL = s2<double>(negposLLStr);
    double negnegLL = s2<double>(negnegLLStr);

    n++;
    if (negnegLL >= negposLL) {
      n_correct++;
    }
    line++;
  }

  FileAccessWrite(out, output);
  double acc = (n_correct/n)*100,
    sn = (n_pos_correct / n_pos)*100,
    sp = (n_correct-n_pos_correct) / (n-n_pos) * 100,
    acc2 = (sn+sp)/2;
  MSG(n << "\t" << n_pos << "\t" << n-n_pos);
  out << std::setw(10) << std::fixed << std::setprecision(2) << acc2
      << std::setw(10) << std::fixed << std::setprecision(2) << sn
      << std::setw(10) << std::fixed << std::setprecision(2) << sp
      << std::setw(10) << std::fixed << std::setprecision(2) << acc
      << ENDL;
}
